multilevel tadpoles
(Ignore probability and estimation for now)
\[ Y_i = b_0 + b_1X_{1i} + ... + \epsilon_i \] This assumes each observation \((i)\) is independent of all other observations. But what if our data were clustered in some way? We might be able to express each observation \(i\) as belonging to a specific cluster, \(j\).
\[ Y_{ij} = b_0 + b_1X_{1ij} + ... + \epsilon_{ij} \]
But what is the problem with this?
We’re starting our unit on multilevel models, which can be thought of as models that “remember” features of clusters of data as they learn about all the clusters. The model will pool information across clusters (e.g., our estimates about cluster A will be informed in part by clusters B, C, and D). This tends to improve estimates about each cluster. Here are some other benefits of multilevel modeling:
Multilevel modeling should be your default approach.
(Still non-Bayesian)
\[\begin{align*} \text{Level 1}&\\ Y_{ij} &= \beta_{0j} + \epsilon_{ij}\\ \text{Level 2}&\\ \beta_{0j} &= \gamma_{00} + U_{0j} \\ U_{0j} &\sim \text{Normal}(0, \tau^2_{00}) \\ \epsilon_{ij} &\sim \text{Normal}(0, \sigma^2) \\ \end{align*}\]
We can rewrite this as the “combined” model:
\[ Y_{ij} = \gamma_{00} + U_{0j} + \epsilon_{ij}\\ \] Level 1 is where you have data that repeats within your grouping or clustering data. Is your cluster classrooms? Then students are Level 1. Is your cluster people? Then observations are Level 1.
Level 2 takes the parameters at Level 1 and decomposes them into a fixed component \((\gamma)\) that reflects the average and, if desired, the individual deviations around that fixed effect \((U)\).
\[\begin{align*} \text{Level 1}&\\ Y_{ij} &= \beta_{0j} + \epsilon_{ij}\\ \text{Level 2}&\\ \beta_{0j} &= \gamma_{00} + U_{0j} \\ U_{0j} &\sim \text{Normal}(0, \tau^2_{00}) \\ \epsilon_{ij} &\sim \text{Normal}(0, \sigma^2) \\ \end{align*}\]
Each person’s score is therefore divided into the “fixed” and “random” components:
\[\begin{align*} Y_{i,16} &= \beta_{0, 16} + \epsilon_{ij} \\ ... \\ B_{0, 16} &= \gamma_{00} + U_{0,16} \end{align*}\]
(Random noise is assumed to vary constantly across observations across participants (in this model).)
iFit study
mean sd 5.5% 94.5% histogram
id 98.61029694 63.7493291 10.0000000 207.000000 ▇▇▇▃▅▅▅▃▃▃▂▂
female 0.57016803 0.4950710 0.0000000 1.000000 ▅▁▁▁▁▁▁▁▁▇
PA.std 0.01438236 1.0241384 -1.6812971 1.751466 ▁▁▃▇▇▃▁▁
day 44.17962096 27.6985612 6.0000000 92.000000 ▇▇▇▇▅▅▅▃▃▃
PA_lag 0.01992587 1.0183833 -1.7204351 1.717036 ▁▂▃▅▇▇▅▃▂▁
NA_lag -0.04575229 0.9833161 -0.8750718 1.990468 ▇▃▂▁▁▁▁▁▁▁▁▁▁▁
steps.pm 0.05424387 0.6298941 -1.0258068 1.011356 ▁▂▇▇▃▁▁▁
steps.pmd 0.02839160 0.7575395 -1.1235951 1.229974 ▁▁▇▇▁▁▁▁▁▁▁▁▁▁
NA.std -0.07545093 0.9495660 -1.0484293 1.811061 ▁▁▁▇▂▁▁▁▁▁set.seed(114)
sample_id = sample(unique(d$id), replace=F, size=6)
d %>%
filter(id %in% sample_id) %>%
mutate(id = factor(id)) %>%
ggplot(aes(x = day, y = PA.std, group = id, fill = id)) +
geom_point(aes(color = id)) +
geom_line( aes(color = id), alpha=.5) +
labs(y="Positive Affect") +
facet_wrap(~id) +
theme(legend.position = "none")What if we wanted to estimate each person’s positive affect? One method would be to simply average scores for each person. This is what we refer to as COMPLETE POOILNG.
Obviously this is not great. People differ, and we probably know we shouldn’t assume each participant has the same level of positive affect. (Although… isn’t this how we talk about people in an independent samples t test?) I think we can acknowledge that we lose a lot of information this way.
Another option would be to treat each person as a group and model scores as a function of group. This is what would be called NO POOLING.
\[\begin{align*} \text{PA.std}_i &\sim \text{Normal}(\mu_i,\sigma) \\ \mu_i &= \alpha_{\text{id}[i]} \\ \alpha_j &\sim \text{Normal}(0, 1.5) \text{ for }j=1,...,239 \\ \sigma &\sim \text{Exponential}(1) \end{align*}\]
Family: gaussian
Links: mu = identity; sigma = identity
Formula: PA.std ~ 0 + a
a ~ 0 + id.f
Data: d (Number of observations: 13033)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
a_id.f1 0.08 0.06 -0.04 0.21 1.00 11270 2736
a_id.f2 1.94 0.06 1.82 2.06 1.00 10974 2899
a_id.f3 -0.03 0.06 -0.15 0.09 1.00 8695 3022
a_id.f4 -0.31 0.06 -0.43 -0.20 1.00 8990 2238
a_id.f5 1.12 0.14 0.84 1.40 1.00 12665 2989
a_id.f6 -0.53 0.07 -0.66 -0.39 1.00 7837 2500
a_id.f7 0.19 0.07 0.05 0.32 1.00 10511 2437
a_id.f8 2.13 0.06 2.00 2.25 1.00 10646 2581
a_id.f9 1.04 0.08 0.89 1.20 1.00 9015 2893
a_id.f10 0.51 0.06 0.39 0.62 1.00 9612 2709
a_id.f11 0.09 0.06 -0.03 0.21 1.00 9234 2953
a_id.f12 0.97 0.06 0.85 1.09 1.00 8591 2734
a_id.f13 0.35 0.06 0.23 0.47 1.00 10446 3255
a_id.f14 -0.12 0.06 -0.24 0.01 1.00 10432 2906
a_id.f15 -0.48 0.06 -0.60 -0.36 1.00 10630 2845
a_id.f16 0.79 0.07 0.67 0.92 1.00 9938 2910
a_id.f17 -0.07 0.07 -0.19 0.06 1.00 9661 2624
a_id.f20 -0.60 0.06 -0.73 -0.48 1.00 9270 2729
a_id.f23 -0.58 0.06 -0.70 -0.46 1.00 10089 3045
a_id.f24 1.97 0.06 1.86 2.09 1.00 7588 2825
a_id.f25 0.29 0.06 0.17 0.41 1.00 8913 2177
a_id.f26 -0.43 0.06 -0.55 -0.31 1.00 9631 2987
a_id.f27 2.02 0.06 1.89 2.14 1.00 8945 2670
a_id.f28 -1.83 0.07 -1.96 -1.71 1.00 9091 2541
a_id.f29 -0.09 0.06 -0.22 0.03 1.00 12304 2658
a_id.f30 0.92 0.06 0.80 1.04 1.00 11099 2580
a_id.f31 -0.48 0.06 -0.60 -0.36 1.00 10211 2559
a_id.f32 1.34 0.06 1.21 1.46 1.00 9726 2749
a_id.f33 -0.78 0.06 -0.90 -0.66 1.00 10981 2697
a_id.f34 0.30 0.06 0.19 0.42 1.00 8492 2520
a_id.f35 -0.15 0.06 -0.27 -0.02 1.00 10766 2704
a_id.f36 0.09 0.06 -0.03 0.21 1.00 13422 2339
a_id.f37 0.73 0.06 0.61 0.85 1.00 11831 2532
a_id.f38 -1.52 0.06 -1.64 -1.41 1.00 8925 2423
a_id.f39 -0.15 0.06 -0.27 -0.04 1.00 8167 2733
a_id.f40 0.25 0.06 0.12 0.37 1.00 10679 2932
a_id.f41 0.78 0.06 0.65 0.90 1.00 8825 2950
a_id.f42 -0.45 0.06 -0.58 -0.33 1.00 9401 2473
a_id.f43 0.43 0.08 0.28 0.58 1.00 9432 2655
a_id.f44 -0.52 0.06 -0.64 -0.39 1.00 10960 2713
a_id.f46 0.18 0.07 0.05 0.31 1.00 10528 2871
a_id.f47 -0.52 0.06 -0.64 -0.41 1.00 10959 2772
a_id.f48 -0.22 0.07 -0.35 -0.10 1.00 7910 2716
a_id.f49 -1.36 0.06 -1.48 -1.24 1.00 8662 2579
a_id.f50 -0.79 0.06 -0.91 -0.68 1.00 9687 2657
a_id.f51 1.77 0.06 1.65 1.90 1.00 8890 2653
a_id.f52 -0.11 0.06 -0.23 0.01 1.00 10873 2560
a_id.f53 0.31 0.06 0.19 0.42 1.00 9033 2564
a_id.f54 0.23 0.06 0.11 0.35 1.00 10638 2764
a_id.f55 1.74 0.06 1.63 1.86 1.00 10136 2617
a_id.f56 -0.21 0.09 -0.37 -0.04 1.00 9053 2606
a_id.f58 -1.44 0.08 -1.60 -1.28 1.00 9842 2312
a_id.f60 0.34 0.06 0.22 0.45 1.00 10447 2865
a_id.f64 0.88 0.07 0.75 1.01 1.00 10141 2394
a_id.f65 0.78 0.07 0.64 0.91 1.00 11076 2488
a_id.f67 -0.08 0.06 -0.21 0.04 1.00 11294 2590
a_id.f68 0.04 0.06 -0.09 0.16 1.00 9932 2665
a_id.f69 -0.81 0.06 -0.93 -0.69 1.00 8488 2789
a_id.f70 0.21 0.06 0.08 0.33 1.00 8768 2299
a_id.f72 -1.21 0.06 -1.32 -1.10 1.00 9975 2808
a_id.f73 -0.12 0.06 -0.24 -0.00 1.00 10602 2661
a_id.f75 0.01 0.06 -0.11 0.12 1.00 9736 2519
a_id.f76 -0.78 0.12 -1.02 -0.53 1.00 9026 2475
a_id.f80 -0.94 0.10 -1.14 -0.75 1.00 10759 2564
a_id.f81 0.29 0.06 0.17 0.41 1.00 10492 2699
a_id.f82 0.48 0.06 0.36 0.59 1.01 8959 2513
a_id.f83 -1.63 0.06 -1.74 -1.51 1.00 10706 2766
a_id.f85 -0.48 0.06 -0.61 -0.36 1.00 8873 2712
a_id.f86 -0.24 0.08 -0.41 -0.08 1.00 9621 2851
a_id.f88 1.10 0.07 0.95 1.24 1.00 11458 2660
a_id.f91 -0.23 0.06 -0.35 -0.12 1.00 12793 2750
a_id.f92 -0.39 0.06 -0.51 -0.27 1.00 7678 2507
a_id.f94 0.02 0.06 -0.10 0.13 1.00 7897 2806
a_id.f95 -0.01 0.07 -0.15 0.13 1.00 11740 2820
a_id.f96 0.34 0.06 0.22 0.46 1.00 8942 2517
a_id.f97 -0.50 0.06 -0.62 -0.38 1.00 11475 2683
a_id.f98 -0.03 0.06 -0.15 0.09 1.00 12370 2673
a_id.f99 -0.46 0.07 -0.60 -0.33 1.00 9260 2712
a_id.f100 -1.44 0.06 -1.56 -1.31 1.00 12439 2672
a_id.f101 -0.67 0.06 -0.79 -0.56 1.00 10542 2914
a_id.f102 -1.31 0.06 -1.43 -1.19 1.00 10056 2485
a_id.f104 1.34 0.06 1.22 1.46 1.00 11217 2417
a_id.f105 -1.10 0.06 -1.22 -0.97 1.00 13252 2519
a_id.f106 0.49 0.06 0.38 0.61 1.00 7821 2608
a_id.f107 0.57 0.06 0.45 0.69 1.00 9225 2433
a_id.f109 1.66 0.06 1.55 1.78 1.00 9725 2928
a_id.f110 -0.05 0.06 -0.17 0.07 1.00 8123 2753
a_id.f111 0.55 0.06 0.44 0.67 1.00 7923 2988
a_id.f112 -0.58 0.15 -0.87 -0.29 1.00 9583 2744
a_id.f113 0.24 0.06 0.12 0.35 1.00 10205 2840
a_id.f114 -0.70 0.07 -0.84 -0.56 1.00 9191 2473
a_id.f115 0.27 0.08 0.10 0.43 1.01 10034 2335
a_id.f116 0.04 0.06 -0.08 0.16 1.00 8985 2855
a_id.f117 0.45 0.12 0.20 0.69 1.00 10768 2864
a_id.f119 0.07 0.12 -0.16 0.29 1.00 8497 2771
a_id.f120 -0.21 0.08 -0.36 -0.05 1.00 9532 2146
a_id.f121 1.22 0.06 1.10 1.34 1.00 12291 2996
a_id.f123 -0.18 0.06 -0.30 -0.06 1.00 10283 2723
a_id.f125 -0.26 0.07 -0.39 -0.13 1.00 9338 2875
a_id.f127 -1.00 0.06 -1.12 -0.88 1.00 11082 2563
a_id.f128 -1.11 0.06 -1.23 -0.98 1.00 10116 2529
a_id.f129 0.57 0.10 0.39 0.75 1.00 8745 3141
a_id.f132 0.10 0.06 -0.02 0.22 1.00 10010 2707
a_id.f133 0.02 0.06 -0.11 0.14 1.00 8711 2583
a_id.f134 0.21 0.07 0.08 0.34 1.00 10237 2613
a_id.f135 1.03 0.21 0.62 1.44 1.00 11717 2849
a_id.f136 -0.18 0.07 -0.31 -0.05 1.00 11568 2788
a_id.f137 1.66 0.08 1.51 1.83 1.01 10422 2345
a_id.f138 -1.68 0.07 -1.81 -1.55 1.00 10819 2759
a_id.f139 0.14 0.07 -0.00 0.27 1.00 9378 2591
a_id.f140 0.06 0.07 -0.07 0.19 1.00 11997 2515
a_id.f142 1.34 0.07 1.20 1.48 1.00 9723 2380
a_id.f143 -1.15 0.07 -1.28 -1.02 1.00 9321 2730
a_id.f144 -0.26 0.12 -0.50 -0.02 1.00 9379 2805
a_id.f145 -1.34 0.08 -1.50 -1.17 1.00 10669 2879
a_id.f147 0.01 0.09 -0.17 0.20 1.00 10415 2557
a_id.f148 -0.05 0.07 -0.19 0.08 1.00 9780 2429
a_id.f149 0.90 0.07 0.76 1.04 1.00 12587 2651
a_id.f150 1.55 0.07 1.42 1.69 1.00 13210 2557
a_id.f151 0.69 0.07 0.55 0.83 1.00 13447 2697
a_id.f152 -0.83 0.07 -0.97 -0.68 1.00 9650 3028
a_id.f153 -1.35 0.07 -1.50 -1.21 1.00 9672 2439
a_id.f154 -0.77 0.08 -0.92 -0.63 1.00 10363 2583
a_id.f155 -0.39 0.08 -0.54 -0.24 1.00 9490 2736
a_id.f157 0.76 0.07 0.62 0.91 1.00 12031 2628
a_id.f158 0.76 0.08 0.61 0.92 1.00 11619 2615
a_id.f160 -0.17 0.08 -0.33 -0.00 1.00 8099 2270
a_id.f162 0.36 0.08 0.20 0.52 1.00 10579 2518
a_id.f163 -0.24 0.07 -0.38 -0.09 1.00 7979 2688
a_id.f164 -1.39 0.07 -1.53 -1.25 1.00 10235 2906
a_id.f165 1.12 0.07 0.98 1.26 1.00 11826 2540
a_id.f166 0.04 0.07 -0.10 0.19 1.00 9628 2716
a_id.f167 -0.46 0.07 -0.61 -0.32 1.00 11116 2846
a_id.f168 -0.26 0.07 -0.40 -0.12 1.00 8970 2591
a_id.f169 0.36 0.08 0.21 0.52 1.01 9863 2360
a_id.f170 -1.38 0.08 -1.53 -1.23 1.00 10634 2948
a_id.f171 -0.52 0.07 -0.67 -0.37 1.00 10290 2607
a_id.f173 -0.01 0.09 -0.18 0.16 1.00 8169 2938
a_id.f174 -0.12 0.09 -0.29 0.04 1.00 12217 2605
a_id.f176 0.43 0.10 0.22 0.63 1.00 8486 2533
a_id.f177 -0.04 0.10 -0.23 0.15 1.01 8429 2566
a_id.f178 -0.19 0.08 -0.34 -0.04 1.00 12379 2182
a_id.f181 0.68 0.08 0.52 0.84 1.00 10100 2429
a_id.f182 0.40 0.08 0.24 0.56 1.00 10363 2484
a_id.f183 0.61 0.08 0.46 0.76 1.00 11405 2813
a_id.f184 -1.61 0.08 -1.77 -1.45 1.00 10245 2776
a_id.f185 0.25 0.08 0.10 0.40 1.00 9114 2973
a_id.f186 -0.62 0.09 -0.80 -0.45 1.00 8718 2540
a_id.f188 -0.20 0.08 -0.36 -0.04 1.00 10469 2604
a_id.f189 0.26 0.10 0.07 0.45 1.00 11539 2702
a_id.f190 -1.72 0.08 -1.89 -1.56 1.00 8338 2445
a_id.f191 -0.06 0.09 -0.22 0.11 1.00 11183 2658
a_id.f194 0.02 0.09 -0.16 0.20 1.00 9761 2193
a_id.f195 1.43 0.08 1.26 1.60 1.00 8897 3002
a_id.f196 0.02 0.10 -0.17 0.20 1.00 7595 2485
a_id.f197 0.64 0.09 0.47 0.81 1.00 11572 2250
a_id.f198 -1.65 0.09 -1.82 -1.46 1.00 8803 3035
a_id.f199 0.46 0.09 0.28 0.63 1.00 9935 2806
a_id.f201 0.68 0.08 0.53 0.82 1.00 8205 2792
a_id.f202 -0.47 0.09 -0.64 -0.29 1.00 10762 3021
a_id.f203 -1.10 0.10 -1.30 -0.90 1.00 10729 2514
a_id.f204 -1.20 0.10 -1.40 -0.99 1.00 9415 2508
a_id.f205 -0.00 0.10 -0.19 0.19 1.00 9123 2452
a_id.f206 -0.18 0.09 -0.35 -0.01 1.00 9946 2997
a_id.f207 -0.41 0.10 -0.60 -0.22 1.00 12955 2891
a_id.f208 0.35 0.10 0.15 0.55 1.00 12311 2147
a_id.f209 -0.47 0.10 -0.67 -0.28 1.00 9681 2707
a_id.f211 -1.27 0.08 -1.43 -1.11 1.00 9841 2794
a_id.f212 0.12 0.10 -0.06 0.32 1.00 10195 2916
a_id.f213 0.17 0.12 -0.07 0.41 1.00 11124 2958
a_id.f214 0.42 0.11 0.21 0.64 1.00 8942 2451
a_id.f216 -0.06 0.13 -0.32 0.19 1.00 11152 2525
a_id.f217 0.03 0.10 -0.17 0.23 1.00 10467 2976
a_id.f219 0.07 0.10 -0.14 0.27 1.00 9617 2229
a_id.f220 1.01 0.11 0.79 1.22 1.00 9917 2500
a_id.f221 -0.26 0.11 -0.48 -0.04 1.00 11306 2733
a_id.f222 0.16 0.10 -0.04 0.37 1.00 8358 2780
a_id.f223 -0.12 0.12 -0.36 0.12 1.00 9828 2595
a_id.f224 -0.24 0.12 -0.48 -0.01 1.00 9876 2741
a_id.f225 -0.09 0.12 -0.32 0.14 1.00 9183 2742
a_id.f226 -0.74 0.14 -1.02 -0.46 1.00 11079 2789
a_id.f227 -0.06 0.15 -0.37 0.25 1.00 8192 2470
a_id.f228 0.44 0.14 0.17 0.72 1.00 9356 2342
a_id.f229 0.52 0.14 0.24 0.80 1.00 10801 2699
a_id.f230 -0.42 0.14 -0.71 -0.14 1.00 10968 2590
a_id.f231 0.59 0.16 0.28 0.90 1.00 10388 2354
a_id.f232 0.60 0.12 0.37 0.83 1.00 10335 2461
a_id.f233 0.39 0.12 0.14 0.63 1.01 9310 2710
a_id.f234 0.90 0.17 0.55 1.23 1.00 9510 2748
a_id.f236 -1.32 0.14 -1.58 -1.06 1.00 9147 2398
a_id.f237 -0.63 0.15 -0.91 -0.34 1.00 9507 2838
a_id.f238 1.04 0.15 0.76 1.32 1.00 10593 2814
a_id.f239 0.55 0.18 0.20 0.90 1.01 10202 2706
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.60 0.00 0.59 0.61 1.00 7117 2734
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).This is inefficient, in that the model treat each person as entirely separate. Let’s try a PARTIAL POOLING model.
\[\begin{align*} \text{PA.std}_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha_{\text{id}[i]} \\ \alpha_j &\sim \text{Normal}(\bar{\alpha}, \sigma_{\alpha}) \text{ for j in 1...239}\\ \bar{\alpha} &\sim \text{Normal}(0, 1.5)\\ \sigma_{\alpha} &\sim \text{Exponential}(1) \\ \sigma &\sim \text{Exponential}(1) \\ \end{align*}\]
m2 <- brm(
data=d,
family=gaussian,
PA.std ~ 1 + (1 | id),
prior = c( prior(normal(0, 1.5), class=Intercept),
prior(exponential(1), class=sd),
prior(exponential(1), class=sigma)),
iter=4000, warmup=1000, chains=4, cores=4, seed=9, #running this longer so it mixes
control = list(adapt_delta =.99),
file=here("files/models/m71.2")
)Despite the additional complexity, we’ll see that our partial pooling model has both a better fit and fewer effective parameters (i.e., is less complex) than our no pooling model.
Let’s visualize the differences in these.
nd1 = distinct(d, id.f)
post1 = epred_draws(m1, nd1)
nd2 = distinct(d, id)
post2 = epred_draws(m2, nd2)
p1 = post1 %>%
ggplot( aes(y=.epred, x=id.f) ) +
stat_gradientinterval() +
scale_x_discrete(labels=NULL, breaks=NULL) +
labs(x="id", y="PA.std", title = "no pooling")
p2 = post2 %>%
mutate(id=as.factor(id)) %>%
ggplot( aes(y=.epred, x=id) ) +
stat_gradientinterval() +
scale_x_discrete(labels=NULL, breaks=NULL) +
labs(x="id", y="PA.std", title = "partial pooling")
p1 / p2
means1 = post1 %>%
mean_qi(.epred)
means2 = post2 %>%
mean_qi(.epred) %>%
mutate(id=as.factor(id))
means1 %>%
ggplot( aes(x=id.f, y=.epred)) +
geom_hline( aes(yintercept=mean(.epred)),
linetype="dashed") +
geom_point( aes(color="no pooling") ) +
geom_point( aes(x=id, color="partial pooling"),
data=means2,
size=2,
alpha=.4) +
scale_color_manual( values=c("#e07a5f", "#1c5253") ) +
scale_x_discrete(breaks=NULL) +
labs(x="id", y="PA.std")+
theme(legend.position = "top")The amount of pooling is affected by the number of data points – lots of data for each cluster = less pooling.
d_small = d %>%
with_groups(id,
sample_n,
size = 3)
m1s <- brm(
data=d_small,
family=gaussian,
bf(PA.std ~ 0 + a,
a ~ 0 + id.f,
nl = TRUE),
prior = c( prior(normal(0, 1.5), class=b, nlpar=a),
prior(exponential(1), class=sigma)),
iter=2000, warmup=1000, chains=4, cores=4, seed=9,
file=here("files/models/m71.1s")
)
m2s <- brm(
data=d_small,
family=gaussian,
PA.std ~ 1 + (1 | id),
prior = c( prior(normal(0, 1.5), class=Intercept),
prior(exponential(1), class=sd),
prior(exponential(1), class=sigma)),
iter=4000, warmup=1000, chains=4, cores=4, seed=9, #running this longer so it mixes
file=here("files/models/m71.2s")
)
nd1 = distinct(d_small, id.f)
post1 = epred_draws(m1s, nd1)
nd2 = distinct(d_small, id)
post2 = epred_draws(m2s, nd2)
means1 = post1 %>%
mean_qi(.epred)
means2 = post2 %>%
mean_qi(.epred) %>%
mutate(id=as.factor(id))
means1 %>%
ggplot( aes(x=id.f, y=.epred)) +
geom_hline( aes(yintercept=mean(.epred)),
linetype="dashed") +
geom_point( aes(color="no pooling") ) +
geom_point( aes(x=id, color="partial pooling"),
data=means2,
size=2,
alpha=.4) +
scale_color_manual( values=c("#e07a5f", "#1c5253") ) +
scale_x_discrete(breaks=NULL) +
labs(x="id", y="PA.std", title="3 observations per person")+
theme(legend.position = "bottom")
Family: gaussian
Links: mu = identity; sigma = identity
Formula: PA.std ~ 1 + (1 | id)
Data: d (Number of observations: 13033)
Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
total post-warmup draws = 12000
Multilevel Hyperparameters:
~id (Number of levels: 193)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.80 0.04 0.73 0.89 1.02 431 563
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.00 0.06 -0.11 0.12 1.05 123 352
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.60 0.00 0.59 0.61 1.00 15924 8655
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Wait, where are all my estimates?
Let’s look closer at these models. How many parameters does each model have?
What additional parameters?
m1 has a unique intercept for each participant and a standard deviation of scores (1 \(\sigma\)).
m2 is estimating all of that plus a grand mean intercept and the variability of means (\(\sigma_M\)).
(what’s the extra one? brms lists the intercept twice. *shrug emoji*)
There’s both an estimate for an overall mean (b_Intercept) as well as the deviations of each person’s mean from that overall (r_id[#, Intercept]).
[1] "b_Intercept" "sd_id__Intercept" "sigma"
[4] "Intercept" "r_id[1,Intercept]" "r_id[2,Intercept]"
[7] "r_id[3,Intercept]" "r_id[4,Intercept]" "r_id[5,Intercept]"
[10] "r_id[6,Intercept]" "r_id[7,Intercept]" "r_id[8,Intercept]"
[13] "r_id[9,Intercept]" "r_id[10,Intercept]" "r_id[11,Intercept]"
[16] "r_id[12,Intercept]" "r_id[13,Intercept]" "r_id[14,Intercept]"
[19] "r_id[15,Intercept]" "r_id[16,Intercept]" "r_id[17,Intercept]"
[22] "r_id[20,Intercept]" "r_id[23,Intercept]" "r_id[24,Intercept]"
[25] "r_id[25,Intercept]" "r_id[26,Intercept]" "r_id[27,Intercept]"
[28] "r_id[28,Intercept]" "r_id[29,Intercept]" "r_id[30,Intercept]"
[31] "r_id[31,Intercept]" "r_id[32,Intercept]" "r_id[33,Intercept]"
[34] "r_id[34,Intercept]" "r_id[35,Intercept]" "r_id[36,Intercept]"
[37] "r_id[37,Intercept]" "r_id[38,Intercept]" "r_id[39,Intercept]"
[40] "r_id[40,Intercept]" "r_id[41,Intercept]" "r_id[42,Intercept]"
[43] "r_id[43,Intercept]" "r_id[44,Intercept]" "r_id[46,Intercept]"
[46] "r_id[47,Intercept]" "r_id[48,Intercept]" "r_id[49,Intercept]"
[49] "r_id[50,Intercept]" "r_id[51,Intercept]" "r_id[52,Intercept]"
[52] "r_id[53,Intercept]" "r_id[54,Intercept]" "r_id[55,Intercept]"
[55] "r_id[56,Intercept]" "r_id[58,Intercept]" "r_id[60,Intercept]"
[58] "r_id[64,Intercept]" "r_id[65,Intercept]" "r_id[67,Intercept]"
[61] "r_id[68,Intercept]" "r_id[69,Intercept]" "r_id[70,Intercept]"
[64] "r_id[72,Intercept]" "r_id[73,Intercept]" "r_id[75,Intercept]"
[67] "r_id[76,Intercept]" "r_id[80,Intercept]" "r_id[81,Intercept]"
[70] "r_id[82,Intercept]" "r_id[83,Intercept]" "r_id[85,Intercept]"
[73] "r_id[86,Intercept]" "r_id[88,Intercept]" "r_id[91,Intercept]"
[76] "r_id[92,Intercept]" "r_id[94,Intercept]" "r_id[95,Intercept]"
[79] "r_id[96,Intercept]" "r_id[97,Intercept]" "r_id[98,Intercept]"
[82] "r_id[99,Intercept]" "r_id[100,Intercept]" "r_id[101,Intercept]"
[85] "r_id[102,Intercept]" "r_id[104,Intercept]" "r_id[105,Intercept]"
[88] "r_id[106,Intercept]" "r_id[107,Intercept]" "r_id[109,Intercept]"
[91] "r_id[110,Intercept]" "r_id[111,Intercept]" "r_id[112,Intercept]"
[94] "r_id[113,Intercept]" "r_id[114,Intercept]" "r_id[115,Intercept]"
[97] "r_id[116,Intercept]" "r_id[117,Intercept]" "r_id[119,Intercept]"
[100] "r_id[120,Intercept]" "r_id[121,Intercept]" "r_id[123,Intercept]"
[103] "r_id[125,Intercept]" "r_id[127,Intercept]" "r_id[128,Intercept]"
[106] "r_id[129,Intercept]" "r_id[132,Intercept]" "r_id[133,Intercept]"
[109] "r_id[134,Intercept]" "r_id[135,Intercept]" "r_id[136,Intercept]"
[112] "r_id[137,Intercept]" "r_id[138,Intercept]" "r_id[139,Intercept]"
[115] "r_id[140,Intercept]" "r_id[142,Intercept]" "r_id[143,Intercept]"
[118] "r_id[144,Intercept]" "r_id[145,Intercept]" "r_id[147,Intercept]"
[121] "r_id[148,Intercept]" "r_id[149,Intercept]" "r_id[150,Intercept]"
[124] "r_id[151,Intercept]" "r_id[152,Intercept]" "r_id[153,Intercept]"
[127] "r_id[154,Intercept]" "r_id[155,Intercept]" "r_id[157,Intercept]"
[130] "r_id[158,Intercept]" "r_id[160,Intercept]" "r_id[162,Intercept]"
[133] "r_id[163,Intercept]" "r_id[164,Intercept]" "r_id[165,Intercept]"
[136] "r_id[166,Intercept]" "r_id[167,Intercept]" "r_id[168,Intercept]"
[139] "r_id[169,Intercept]" "r_id[170,Intercept]" "r_id[171,Intercept]"
[142] "r_id[173,Intercept]" "r_id[174,Intercept]" "r_id[176,Intercept]"
[145] "r_id[177,Intercept]" "r_id[178,Intercept]" "r_id[181,Intercept]"
[148] "r_id[182,Intercept]" "r_id[183,Intercept]" "r_id[184,Intercept]"
[151] "r_id[185,Intercept]" "r_id[186,Intercept]" "r_id[188,Intercept]"
[154] "r_id[189,Intercept]" "r_id[190,Intercept]" "r_id[191,Intercept]"
[157] "r_id[194,Intercept]" "r_id[195,Intercept]" "r_id[196,Intercept]"
[160] "r_id[197,Intercept]" "r_id[198,Intercept]" "r_id[199,Intercept]"
[163] "r_id[201,Intercept]" "r_id[202,Intercept]" "r_id[203,Intercept]"
[166] "r_id[204,Intercept]" "r_id[205,Intercept]" "r_id[206,Intercept]"
[169] "r_id[207,Intercept]" "r_id[208,Intercept]" "r_id[209,Intercept]"
[172] "r_id[211,Intercept]" "r_id[212,Intercept]" "r_id[213,Intercept]"
[175] "r_id[214,Intercept]" "r_id[216,Intercept]" "r_id[217,Intercept]"
[178] "r_id[219,Intercept]" "r_id[220,Intercept]" "r_id[221,Intercept]"
[181] "r_id[222,Intercept]" "r_id[223,Intercept]" "r_id[224,Intercept]"
[184] "r_id[225,Intercept]" "r_id[226,Intercept]" "r_id[227,Intercept]"
[187] "r_id[228,Intercept]" "r_id[229,Intercept]" "r_id[230,Intercept]"
[190] "r_id[231,Intercept]" "r_id[232,Intercept]" "r_id[233,Intercept]"
[193] "r_id[234,Intercept]" "r_id[236,Intercept]" "r_id[237,Intercept]"
[196] "r_id[238,Intercept]" "r_id[239,Intercept]" "lprior"
[199] "lp__" "accept_stat__" "stepsize__"
[202] "treedepth__" "n_leapfrog__" "divergent__"
[205] "energy__" You can combine these to get each participant’s model-implied mean.
# A tibble: 2,316,000 × 8
# Groups: id, Intercept [193]
.chain .iteration .draw b_Intercept id Intercept r_id PA
<int> <int> <int> <dbl> <int> <chr> <dbl> <dbl>
1 1 1 1 0.0269 1 Intercept 0.0819 0.109
2 1 1 1 0.0269 2 Intercept 1.88 1.91
3 1 1 1 0.0269 3 Intercept -0.188 -0.161
4 1 1 1 0.0269 4 Intercept -0.380 -0.353
5 1 1 1 0.0269 5 Intercept 1.17 1.20
6 1 1 1 0.0269 6 Intercept -0.433 -0.406
7 1 1 1 0.0269 7 Intercept 0.102 0.129
8 1 1 1 0.0269 8 Intercept 2.10 2.12
9 1 1 1 0.0269 9 Intercept 0.889 0.916
10 1 1 1 0.0269 10 Intercept 0.339 0.366
# ℹ 2,315,990 more rowsOr (when in doubt) just use the expected predictions function (see how they’re the same!):
# A tibble: 2,316,000 × 6
# Groups: id, .row [193]
id .row .chain .iteration .draw .epred
<int> <int> <int> <int> <int> <dbl>
1 1 1 NA NA 1 0.109
2 1 1 NA NA 2 0.208
3 1 1 NA NA 3 -0.0119
4 1 1 NA NA 4 0.203
5 1 1 NA NA 5 -0.0400
6 1 1 NA NA 6 0.132
7 1 1 NA NA 7 0.0361
8 1 1 NA NA 8 0.116
9 1 1 NA NA 9 0.0561
10 1 1 NA NA 10 0.0557
# ℹ 2,315,990 more rowsLet’s see too how this model has decomposed the within- and between-person variance.
sd var
1.024 1.049 # A tibble: 12,000 × 5
.chain .iteration .draw sd_id__Intercept sigma
<int> <int> <int> <dbl> <dbl>
1 1 1 1 0.886 0.600
2 1 2 2 0.872 0.604
3 1 3 3 0.871 0.600
4 1 4 4 0.887 0.602
5 1 5 5 0.892 0.601
6 1 6 6 0.886 0.606
7 1 7 7 0.846 0.599
8 1 8 8 0.847 0.601
9 1 9 9 0.849 0.605
10 1 10 10 0.844 0.607
# ℹ 11,990 more rows\[\begin{align*} \text{PA.std}_{ij} &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha_{id[i]} \\ \alpha_{id[i]} &\sim \text{Normal}(\bar{\alpha}, 1) \text{ for i in } 1...239\\ \bar{\alpha} &\sim \text{Normal}(.50,.25) \\ \sigma &\sim \text{Exponential}(1) \\ \end{align*}\]
Another way to write out the expected value of \(\mu_j\) is:
\[ \mu_j = \bar{\alpha} + \alpha_{i[ID]} \]
Where \(\bar{\alpha}\) is the mean of the means and \(\alpha_{i[ID]}\) is the deviation of each person from the grand mean. There will be one \(\bar{\alpha}\) for the entire sample, but a different \(\alpha_{i[ID]}\) for each person. Because there are multiple values for each person, we can calculate the variability of those \(\alpha\)’s, as well as the residual variability within each person.
\[\begin{align*} \text{PA.std}_{ij} &\sim \text{Normal}(\mu_j, \sigma) \\ \mu_i &= \bar{\alpha} + \alpha_{id[j]} \\ \bar{\alpha} &\sim \text{Normal}(.50,.25) \\ \sigma &\sim \text{Exponential}(1) \\ \sigma_{\alpha} &\sim \text{Exponential}(1) \\ \end{align*}\]
In this reparameterization, we no longer need a prior for each of the person-means. That’s because they must have a mean of 0. Instead, we only need to estimate the variability \((\sigma_{\alpha})\) of these deviations.
In other words, this model is analogous to an ANOVA, in which we calculate both within- and between-group (aka person) variability. Therefore, we’ll have a standard deviation for both levels.
In this figure, we can visualize three approaches to pooling:
actual_means = d %>%
with_groups(id, summarise, p_mean = mean(PA.std)) %>%
mutate()
spread_draws(m2, b_Intercept, r_id[id, term]) %>%
mutate(p_mean = b_Intercept + r_id) %>%
mean_qi(p_mean) %>%
filter(id %in% sample_id) %>%
ggplot( aes(x=id, y=p_mean) ) +
geom_hline( aes(yintercept = mean(d$PA.std)),
linetype="dashed") +
geom_point( size=2, color="#e07a5f") +
geom_point( data=filter(actual_means,
id %in% sample_id)) +
labs(x="id",y="PA.std")
It’s common to write out mulitlevel models using formulas for the different levels. Level 1 is the level of your outcome, or the thing that repeats. Level 2 is the level of your groups.
\[\begin{align*} \text{Level 1} &\\ \text{PA.std}_{ij} &= \alpha_j + \epsilon_{ij} \\ \text{Level 2} &\\ \alpha_j &= \gamma_0 + U_j \end{align*}\]
Some refer to \(U_j\) as a “random” or “varying” effect because it varies across groups. \(\gamma_0\) is therefore a “fixed” or “non-varying” effect because it applies to all groups.
We’ve already seen how we can use spread_draws to get our parameter estimates.
# A tibble: 6 × 6
.chain .iteration .draw b_Intercept sigma sd_id__Intercept
<int> <int> <int> <dbl> <dbl> <dbl>
1 1 1 1 0.0269 0.600 0.886
2 1 2 2 0.0313 0.604 0.872
3 1 3 3 0.0264 0.600 0.871
4 1 4 4 0.0217 0.602 0.887
5 1 5 5 0.0318 0.601 0.892
6 1 6 6 0.0273 0.606 0.886# A tibble: 6 × 6
# Groups: id, term [1]
id term r_id .chain .iteration .draw
<int> <chr> <dbl> <int> <int> <int>
1 1 Intercept 0.0819 1 1 1
2 1 Intercept 0.176 1 2 2
3 1 Intercept -0.0383 1 3 3
4 1 Intercept 0.182 1 4 4
5 1 Intercept -0.0718 1 5 5
6 1 Intercept 0.104 1 6 6We can get expected values (means) for individuals in our sample.
set.seed(9)
sample_id = sample(unique(d$id), replace=F, size=3)
nd = data.frame(id=sample_id, id.f = sample_id)
add_epred_draws(newdata=nd, object=m2) %>%
mean_qi()# A tibble: 3 × 9
id id.f .row .epred .lower .upper .width .point .interval
<int> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 6 6 3 -0.525 -0.664 -0.387 0.95 mean qi
2 60 60 2 0.334 0.213 0.454 0.95 mean qi
3 232 232 1 0.594 0.366 0.820 0.95 mean qi We can get predicted values (scores) for individuals in our sample.
# A tibble: 3 × 9
id id.f .row .prediction .lower .upper .width .point .interval
<int> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 6 6 3 -0.531 -1.72 0.632 0.95 mean qi
2 60 60 2 0.324 -0.886 1.52 0.95 mean qi
3 232 232 1 0.586 -0.600 1.77 0.95 mean qi expected = add_epred_draws(newdata=nd, object=m2)
predicted = add_predicted_draws(newdata=nd, object=m2)
full_join(expected, predicted) %>%
filter(.draw <= 100) %>%
ungroup() %>%
ggplot( aes(x=id, y=.epred)) +
stat_halfeye( aes(fill="expected"),
alpha=.3) +
stat_halfeye( aes(y =.prediction, fill="predicted"),
alpha=.3 ) +
scale_fill_manual( values=c("#e07a5f","#1c5253") ) +
labs(x="id", y="PA.std", fill="draw") +
scale_x_continuous(breaks=sample_id) +
theme(legend.position = "top")Finally, we can get expected values for new individuals
Let’s add time to our model. Because time varies (can change) within person from assessment to assessment, this is a Level 1 variable. Note that I have NOT added a varying component to Level 2 – in other words, I’m stating that the effect of time on positive affect is fixed or identical across participants.
\[\begin{align*} \text{Level 1} &\\ \text{PA.std}_{ij} &= \alpha_i + \beta_i(\text{day}_{ij}) + \epsilon_{ij} \\ \text{Level 2} &\\ \alpha_j &= \gamma_0 + U_{0j} \\ \beta_j &= \gamma_1 \\ \end{align*}\]
m3 <- brm(
data=d,
family=gaussian,
PA.std ~ 1 + day + (1 | id),
prior = c( prior(normal(.50, .25), class=Intercept),
prior(normal(0, 1), class=b), #new prior for new term
prior(exponential(1), class=sd),
prior(exponential(1), class=sigma)),
iter=2000, warmup=1000, chains=4, cores=4, seed=9,
file=here("files/models/m71.3")
) Family: gaussian
Links: mu = identity; sigma = identity
Formula: PA.std ~ 1 + day + (1 | id)
Data: d (Number of observations: 13033)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Multilevel Hyperparameters:
~id (Number of levels: 193)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.80 0.04 0.73 0.88 1.04 110 285
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.13 0.05 0.02 0.23 1.08 47 188
day -0.00 0.00 -0.00 -0.00 1.00 4550 3272
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.60 0.00 0.59 0.60 1.00 6598 2554
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).There are different intercepts for each participant, but not different slopes.
[1] "b_Intercept" "b_day" "sd_id__Intercept"
[4] "sigma" "Intercept" "r_id[1,Intercept]"
[7] "r_id[2,Intercept]" "r_id[3,Intercept]" "r_id[4,Intercept]"
[10] "r_id[5,Intercept]" "r_id[6,Intercept]" "r_id[7,Intercept]"
[13] "r_id[8,Intercept]" "r_id[9,Intercept]" "r_id[10,Intercept]"
[16] "r_id[11,Intercept]" "r_id[12,Intercept]" "r_id[13,Intercept]"
[19] "r_id[14,Intercept]" "r_id[15,Intercept]" "r_id[16,Intercept]"
[22] "r_id[17,Intercept]" "r_id[20,Intercept]" "r_id[23,Intercept]"
[25] "r_id[24,Intercept]" "r_id[25,Intercept]" "r_id[26,Intercept]"
[28] "r_id[27,Intercept]" "r_id[28,Intercept]" "r_id[29,Intercept]"
[31] "r_id[30,Intercept]" "r_id[31,Intercept]" "r_id[32,Intercept]"
[34] "r_id[33,Intercept]" "r_id[34,Intercept]" "r_id[35,Intercept]"
[37] "r_id[36,Intercept]" "r_id[37,Intercept]" "r_id[38,Intercept]"
[40] "r_id[39,Intercept]" "r_id[40,Intercept]" "r_id[41,Intercept]"
[43] "r_id[42,Intercept]" "r_id[43,Intercept]" "r_id[44,Intercept]"
[46] "r_id[46,Intercept]" "r_id[47,Intercept]" "r_id[48,Intercept]"
[49] "r_id[49,Intercept]" "r_id[50,Intercept]" "r_id[51,Intercept]"
[52] "r_id[52,Intercept]" "r_id[53,Intercept]" "r_id[54,Intercept]"
[55] "r_id[55,Intercept]" "r_id[56,Intercept]" "r_id[58,Intercept]"
[58] "r_id[60,Intercept]" "r_id[64,Intercept]" "r_id[65,Intercept]"
[61] "r_id[67,Intercept]" "r_id[68,Intercept]" "r_id[69,Intercept]"
[64] "r_id[70,Intercept]" "r_id[72,Intercept]" "r_id[73,Intercept]"
[67] "r_id[75,Intercept]" "r_id[76,Intercept]" "r_id[80,Intercept]"
[70] "r_id[81,Intercept]" "r_id[82,Intercept]" "r_id[83,Intercept]"
[73] "r_id[85,Intercept]" "r_id[86,Intercept]" "r_id[88,Intercept]"
[76] "r_id[91,Intercept]" "r_id[92,Intercept]" "r_id[94,Intercept]"
[79] "r_id[95,Intercept]" "r_id[96,Intercept]" "r_id[97,Intercept]"
[82] "r_id[98,Intercept]" "r_id[99,Intercept]" "r_id[100,Intercept]"
[85] "r_id[101,Intercept]" "r_id[102,Intercept]" "r_id[104,Intercept]"
[88] "r_id[105,Intercept]" "r_id[106,Intercept]" "r_id[107,Intercept]"
[91] "r_id[109,Intercept]" "r_id[110,Intercept]" "r_id[111,Intercept]"
[94] "r_id[112,Intercept]" "r_id[113,Intercept]" "r_id[114,Intercept]"
[97] "r_id[115,Intercept]" "r_id[116,Intercept]" "r_id[117,Intercept]"
[100] "r_id[119,Intercept]" "r_id[120,Intercept]" "r_id[121,Intercept]"
[103] "r_id[123,Intercept]" "r_id[125,Intercept]" "r_id[127,Intercept]"
[106] "r_id[128,Intercept]" "r_id[129,Intercept]" "r_id[132,Intercept]"
[109] "r_id[133,Intercept]" "r_id[134,Intercept]" "r_id[135,Intercept]"
[112] "r_id[136,Intercept]" "r_id[137,Intercept]" "r_id[138,Intercept]"
[115] "r_id[139,Intercept]" "r_id[140,Intercept]" "r_id[142,Intercept]"
[118] "r_id[143,Intercept]" "r_id[144,Intercept]" "r_id[145,Intercept]"
[121] "r_id[147,Intercept]" "r_id[148,Intercept]" "r_id[149,Intercept]"
[124] "r_id[150,Intercept]" "r_id[151,Intercept]" "r_id[152,Intercept]"
[127] "r_id[153,Intercept]" "r_id[154,Intercept]" "r_id[155,Intercept]"
[130] "r_id[157,Intercept]" "r_id[158,Intercept]" "r_id[160,Intercept]"
[133] "r_id[162,Intercept]" "r_id[163,Intercept]" "r_id[164,Intercept]"
[136] "r_id[165,Intercept]" "r_id[166,Intercept]" "r_id[167,Intercept]"
[139] "r_id[168,Intercept]" "r_id[169,Intercept]" "r_id[170,Intercept]"
[142] "r_id[171,Intercept]" "r_id[173,Intercept]" "r_id[174,Intercept]"
[145] "r_id[176,Intercept]" "r_id[177,Intercept]" "r_id[178,Intercept]"
[148] "r_id[181,Intercept]" "r_id[182,Intercept]" "r_id[183,Intercept]"
[151] "r_id[184,Intercept]" "r_id[185,Intercept]" "r_id[186,Intercept]"
[154] "r_id[188,Intercept]" "r_id[189,Intercept]" "r_id[190,Intercept]"
[157] "r_id[191,Intercept]" "r_id[194,Intercept]" "r_id[195,Intercept]"
[160] "r_id[196,Intercept]" "r_id[197,Intercept]" "r_id[198,Intercept]"
[163] "r_id[199,Intercept]" "r_id[201,Intercept]" "r_id[202,Intercept]"
[166] "r_id[203,Intercept]" "r_id[204,Intercept]" "r_id[205,Intercept]"
[169] "r_id[206,Intercept]" "r_id[207,Intercept]" "r_id[208,Intercept]"
[172] "r_id[209,Intercept]" "r_id[211,Intercept]" "r_id[212,Intercept]"
[175] "r_id[213,Intercept]" "r_id[214,Intercept]" "r_id[216,Intercept]"
[178] "r_id[217,Intercept]" "r_id[219,Intercept]" "r_id[220,Intercept]"
[181] "r_id[221,Intercept]" "r_id[222,Intercept]" "r_id[223,Intercept]"
[184] "r_id[224,Intercept]" "r_id[225,Intercept]" "r_id[226,Intercept]"
[187] "r_id[227,Intercept]" "r_id[228,Intercept]" "r_id[229,Intercept]"
[190] "r_id[230,Intercept]" "r_id[231,Intercept]" "r_id[232,Intercept]"
[193] "r_id[233,Intercept]" "r_id[234,Intercept]" "r_id[236,Intercept]"
[196] "r_id[237,Intercept]" "r_id[238,Intercept]" "r_id[239,Intercept]"
[199] "lprior" "lp__" "accept_stat__"
[202] "stepsize__" "treedepth__" "n_leapfrog__"
[205] "divergent__" "energy__" We can also add covariates at Level 2, or the person-level in this case. We have a binary variable called female that refers to the gender of each person1.
\[\begin{align*} \text{Level 1} &\\ \text{PA.std}_{ij} &= \alpha_i + \beta_i(\text{day}_{ij}) + \epsilon_{ij} \\ \text{Level 2} &\\ \alpha_j &= \gamma_{0\text{female}[j]} + U_{0j} \\ \beta_j &= \gamma_{1} \\ \end{align*}\]
Family: gaussian
Links: mu = identity; sigma = identity
Formula: PA.std ~ 0 + day + female + (1 | id)
Data: d (Number of observations: 13033)
Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
total post-warmup draws = 16000
Multilevel Hyperparameters:
~id (Number of levels: 193)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.77 0.04 0.70 0.85 1.00 632 1163
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
day -0.00 0.00 -0.00 -0.00 1.00 13839 13412
female 0.35 0.07 0.20 0.49 1.01 371 731
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.60 0.00 0.59 0.60 1.00 36313 10598
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1). [1] "b_day" "b_female" "sd_id__Intercept"
[4] "sigma" "r_id[1,Intercept]" "r_id[2,Intercept]"
[7] "r_id[3,Intercept]" "r_id[4,Intercept]" "r_id[5,Intercept]"
[10] "r_id[6,Intercept]" "r_id[7,Intercept]" "r_id[8,Intercept]"
[13] "r_id[9,Intercept]" "r_id[10,Intercept]" "r_id[11,Intercept]"
[16] "r_id[12,Intercept]" "r_id[13,Intercept]" "r_id[14,Intercept]"
[19] "r_id[15,Intercept]" "r_id[16,Intercept]" "r_id[17,Intercept]"
[22] "r_id[20,Intercept]" "r_id[23,Intercept]" "r_id[24,Intercept]"
[25] "r_id[25,Intercept]" "r_id[26,Intercept]" "r_id[27,Intercept]"
[28] "r_id[28,Intercept]" "r_id[29,Intercept]" "r_id[30,Intercept]"
[31] "r_id[31,Intercept]" "r_id[32,Intercept]" "r_id[33,Intercept]"
[34] "r_id[34,Intercept]" "r_id[35,Intercept]" "r_id[36,Intercept]"
[37] "r_id[37,Intercept]" "r_id[38,Intercept]" "r_id[39,Intercept]"
[40] "r_id[40,Intercept]" "r_id[41,Intercept]" "r_id[42,Intercept]"
[43] "r_id[43,Intercept]" "r_id[44,Intercept]" "r_id[46,Intercept]"
[46] "r_id[47,Intercept]" "r_id[48,Intercept]" "r_id[49,Intercept]"
[49] "r_id[50,Intercept]" "r_id[51,Intercept]" "r_id[52,Intercept]"
[52] "r_id[53,Intercept]" "r_id[54,Intercept]" "r_id[55,Intercept]"
[55] "r_id[56,Intercept]" "r_id[58,Intercept]" "r_id[60,Intercept]"
[58] "r_id[64,Intercept]" "r_id[65,Intercept]" "r_id[67,Intercept]"
[61] "r_id[68,Intercept]" "r_id[69,Intercept]" "r_id[70,Intercept]"
[64] "r_id[72,Intercept]" "r_id[73,Intercept]" "r_id[75,Intercept]"
[67] "r_id[76,Intercept]" "r_id[80,Intercept]" "r_id[81,Intercept]"
[70] "r_id[82,Intercept]" "r_id[83,Intercept]" "r_id[85,Intercept]"
[73] "r_id[86,Intercept]" "r_id[88,Intercept]" "r_id[91,Intercept]"
[76] "r_id[92,Intercept]" "r_id[94,Intercept]" "r_id[95,Intercept]"
[79] "r_id[96,Intercept]" "r_id[97,Intercept]" "r_id[98,Intercept]"
[82] "r_id[99,Intercept]" "r_id[100,Intercept]" "r_id[101,Intercept]"
[85] "r_id[102,Intercept]" "r_id[104,Intercept]" "r_id[105,Intercept]"
[88] "r_id[106,Intercept]" "r_id[107,Intercept]" "r_id[109,Intercept]"
[91] "r_id[110,Intercept]" "r_id[111,Intercept]" "r_id[112,Intercept]"
[94] "r_id[113,Intercept]" "r_id[114,Intercept]" "r_id[115,Intercept]"
[97] "r_id[116,Intercept]" "r_id[117,Intercept]" "r_id[119,Intercept]"
[100] "r_id[120,Intercept]" "r_id[121,Intercept]" "r_id[123,Intercept]"
[103] "r_id[125,Intercept]" "r_id[127,Intercept]" "r_id[128,Intercept]"
[106] "r_id[129,Intercept]" "r_id[132,Intercept]" "r_id[133,Intercept]"
[109] "r_id[134,Intercept]" "r_id[135,Intercept]" "r_id[136,Intercept]"
[112] "r_id[137,Intercept]" "r_id[138,Intercept]" "r_id[139,Intercept]"
[115] "r_id[140,Intercept]" "r_id[142,Intercept]" "r_id[143,Intercept]"
[118] "r_id[144,Intercept]" "r_id[145,Intercept]" "r_id[147,Intercept]"
[121] "r_id[148,Intercept]" "r_id[149,Intercept]" "r_id[150,Intercept]"
[124] "r_id[151,Intercept]" "r_id[152,Intercept]" "r_id[153,Intercept]"
[127] "r_id[154,Intercept]" "r_id[155,Intercept]" "r_id[157,Intercept]"
[130] "r_id[158,Intercept]" "r_id[160,Intercept]" "r_id[162,Intercept]"
[133] "r_id[163,Intercept]" "r_id[164,Intercept]" "r_id[165,Intercept]"
[136] "r_id[166,Intercept]" "r_id[167,Intercept]" "r_id[168,Intercept]"
[139] "r_id[169,Intercept]" "r_id[170,Intercept]" "r_id[171,Intercept]"
[142] "r_id[173,Intercept]" "r_id[174,Intercept]" "r_id[176,Intercept]"
[145] "r_id[177,Intercept]" "r_id[178,Intercept]" "r_id[181,Intercept]"
[148] "r_id[182,Intercept]" "r_id[183,Intercept]" "r_id[184,Intercept]"
[151] "r_id[185,Intercept]" "r_id[186,Intercept]" "r_id[188,Intercept]"
[154] "r_id[189,Intercept]" "r_id[190,Intercept]" "r_id[191,Intercept]"
[157] "r_id[194,Intercept]" "r_id[195,Intercept]" "r_id[196,Intercept]"
[160] "r_id[197,Intercept]" "r_id[198,Intercept]" "r_id[199,Intercept]"
[163] "r_id[201,Intercept]" "r_id[202,Intercept]" "r_id[203,Intercept]"
[166] "r_id[204,Intercept]" "r_id[205,Intercept]" "r_id[206,Intercept]"
[169] "r_id[207,Intercept]" "r_id[208,Intercept]" "r_id[209,Intercept]"
[172] "r_id[211,Intercept]" "r_id[212,Intercept]" "r_id[213,Intercept]"
[175] "r_id[214,Intercept]" "r_id[216,Intercept]" "r_id[217,Intercept]"
[178] "r_id[219,Intercept]" "r_id[220,Intercept]" "r_id[221,Intercept]"
[181] "r_id[222,Intercept]" "r_id[223,Intercept]" "r_id[224,Intercept]"
[184] "r_id[225,Intercept]" "r_id[226,Intercept]" "r_id[227,Intercept]"
[187] "r_id[228,Intercept]" "r_id[229,Intercept]" "r_id[230,Intercept]"
[190] "r_id[231,Intercept]" "r_id[232,Intercept]" "r_id[233,Intercept]"
[193] "r_id[234,Intercept]" "r_id[236,Intercept]" "r_id[237,Intercept]"
[196] "r_id[238,Intercept]" "r_id[239,Intercept]" "lprior"
[199] "lp__" "accept_stat__" "stepsize__"
[202] "treedepth__" "n_leapfrog__" "divergent__"
[205] "energy__" set.seed(9)
sample_id = sample(unique(d$id), replace=F, size = 20)
distinct(d, id, female, day) %>%
filter(id %in% sample_id) %>%
add_epred_draws(m4) %>%
mean_qi() %>%
ggplot( aes(x=day, y=.epred, group = as.factor(id), color=as.factor(female))) +
geom_line() +
scale_color_manual(values=c("#1c5253" , "#e07a5f")) +
labs(x="day",y="PA.std", color = "female") +
facet_wrap(~female) +
theme(legend.position = "top")So far, we’ve only added covariates to the estimates of the intercept. But we can have covariates on our slope parameters. These are interactions (the effect of variable X on Y changes as a function of variable Z.) And they introduce additional parameters into our model.
\[\begin{align*} \text{Level 1} &\\ \text{PA.std}_{ij} &= \alpha_i + \beta_i(\text{day}_{ij}) + \epsilon_{ij} \\ \text{Level 2} &\\ \alpha_j &= \gamma_{00} + \gamma_{00}\text{female}_j + U_{0j} \\ \beta_j &= \gamma_{10} + \gamma_{10}\text{female}_j + U_{1j} \\ \end{align*}\]
m5 <- brm(
data=d,
family=gaussian,
PA.std ~ 1 + day*female + (1 + day | id),
prior = c( prior(normal(0, 1), class=Intercept),
prior(normal(0, 1), class=b),
prior(exponential(1), class=sd),
prior(exponential(1), class=sigma)),
iter=5000, warmup=1000, chains=4, cores=4, seed=9,
file=here("files/models/m71.5")
)Note the hyperparameter, cor(Intercept, day)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: PA.std ~ 1 + day * female + (1 + day | id)
Data: d (Number of observations: 13033)
Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
total post-warmup draws = 16000
Multilevel Hyperparameters:
~id (Number of levels: 193)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.79 0.04 0.72 0.88 1.00 2408 5140
sd(day) 0.01 0.00 0.01 0.01 1.00 4026 6741
cor(Intercept,day) -0.29 0.07 -0.43 -0.15 1.00 2667 5616
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -0.16 0.09 -0.32 0.02 1.00 1290 2757
day -0.00 0.00 -0.01 -0.00 1.00 2518 4877
female 0.49 0.11 0.27 0.71 1.00 1228 2512
day:female 0.00 0.00 -0.00 0.00 1.00 2271 5046
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.56 0.00 0.56 0.57 1.00 35997 10764
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).